function mysimulation(fs,t_initial,n_sample)
h = 1/fs;                   % integration steps
t_final = t_initial + h*n_sample
t=t_initial:h:t_final;		% setup time from 0 to t_final
tcnt = t_initial:1/10000:t_final;   % times continuous vector
f=50;						% input frequency
f1=100;                     % input frequency
f2=200;
f3=300;
f4=400;
f5=500;
f6=600;
f7=700;
f8=800;
f9=900;
f10=1000;


n_int = length(t);          % required integration
% dydt(1)=2*pi*f*cos(2*pi*f*t(1)); % Initial acceleration value

input = 2*pi*f*cos(2*pi*f*t)+ 2*pi*f1*cos(2*pi*f1*t) + 2*pi*f2*cos(2*pi*f2*t)...
    + 2*pi*f3*cos(2*pi*f3*t) + 2*pi*f4*cos(2*pi*f4*t);

inputcnt = 2*pi*f*cos(2*pi*f*tcnt)+ 2*pi*f1*cos(2*pi*f1*tcnt) + 2*pi*f2*cos(2*pi*f2*tcnt)...
    + 2*pi*f3*cos(2*pi*f3*tcnt) + 2*pi*f4*cos(2*pi*f4*tcnt);

% input = 2*pi*f*cos(2*pi*f*t)+ 2*pi*f1*cos(2*pi*f1*t) + 2*pi*f2*cos(2*pi*f2*t)...
%     + 2*pi*f3*cos(2*pi*f3*t) + 2*pi*f4*cos(2*pi*f4*t) + 2*pi*f5*cos(2*pi*f5*t)...
%     + 2*pi*f6*cos(2*pi*f6*t) + 2*pi*f7*cos(2*pi*f7*t) + 2*pi*f8*cos(2*pi*f8*t)...
%     + 2*pi*f9*cos(2*pi*f9*t) + 2*pi*f10*cos(2*pi*f10*t); 
figure;
% subplot(2,1,1);
plot(tcnt,inputcnt);
hold on
title('Signal Input');
xlabel('time (seconds)');
% subplot(2,1,2);
stem(t,input);
% title('Sample Input');
% xlabel('time (seconds)');
hold off

syms x
% fx = 2*pi*f*cos(2*pi*f*x) ;
fx = 2*pi*f*cos(2*pi*f*x) + 2*pi*f1*cos(2*pi*f1*x) + 2*pi*f2*cos(2*pi*f2*x)...
    + 2*pi*f3*cos(2*pi*f3*x) + 2*pi*f4*cos(2*pi*f4*x);

% fx = 2*pi*f*cos(2*pi*f*x) + 2*pi*f1*cos(2*pi*f1*x) + 2*pi*f2*cos(2*pi*f2*x)...
%     + 2*pi*f3*cos(2*pi*f3*x) + 2*pi*f4*cos(2*pi*f4*x) + 2*pi*f5*cos(2*pi*f5*x)...
%     + 2*pi*f6*cos(2*pi*f6*x) + 2*pi*f7*cos(2*pi*f7*x) + 2*pi*f8*cos(2*pi*f8*x)...
%     + 2*pi*f9*cos(2*pi*f9*x) + 2*pi*f10*cos(2*pi*f10*x);

%pre alocating value
expected_value = zeros(1,n_int);
y = zeros(1,n_int);

% Initial value
y(1) = 0;                   
expected_value(1) = double(int(fx,x,0,t(1)));

dydt(1)=2*pi*f*cos(2*pi*f*t(1))+ 2*pi*f1*cos(2*pi*f1*t(1)) + 2*pi*f2*cos(2*pi*f2*t(1))...
    + 2*pi*f3*cos(2*pi*f3*t(1)) + 2*pi*f4*cos(2*pi*f4*t(1));

% dydt(1)=2*pi*f*cos(2*pi*f*t(1))+ 2*pi*f1*cos(2*pi*f1*t(1)) + 2*pi*f2*cos(2*pi*f2*t(1))...
%     + 2*pi*f3*cos(2*pi*f3*t(1)) + 2*pi*f4*cos(2*pi*f4*t(1)) + 2*pi*f5*cos(2*pi*f5*t(1))...
%     + 2*pi*f6*cos(2*pi*f6*t(1)) + 2*pi*f7*cos(2*pi*f7*t(1)) + 2*pi*f8*cos(2*pi*f8*t(1))...
%     + 2*pi*f9*cos(2*pi*f9*t(1)) + 2*pi*f10*cos(2*pi*f10*t(1));   % Initial acceleration value

for i=2:n_int
    %compute integration
%     dydt(2)=2*pi*f*cos(2*pi*f*t(i));

    dydt(2)=2*pi*f*cos(2*pi*f*t(i)) + 2*pi*f1*cos(2*pi*f1*t(i)) + 2*pi*f2*cos(2*pi*f2*t(i))...
    + 2*pi*f3*cos(2*pi*f3*t(i)) + 2*pi*f4*cos(2*pi*f4*t(i));

%     dydt(2)=2*pi*f*cos(2*pi*f*t(i)) + 2*pi*f1*cos(2*pi*f1*t(i)) + 2*pi*f2*cos(2*pi*f2*t(i))...
%     + 2*pi*f3*cos(2*pi*f3*t(i)) + 2*pi*f4*cos(2*pi*f4*t(i)) + 2*pi*f5*cos(2*pi*f5*t(i))...
%     + 2*pi*f6*cos(2*pi*f6*t(i)) + 2*pi*f7*cos(2*pi*f7*t(i)) + 2*pi*f8*cos(2*pi*f8*t(i))...
%     + 2*pi*f9*cos(2*pi*f9*t(i)) + 2*pi*f10*cos(2*pi*f10*t(i));   % Input signal

    y(i) = y(i-1) + (h/2)*(dydt(2)+dydt(1));
    %save value to next iteration
    dydt(1)=dydt(2);
    
    expected_value(i) = double(int(fx,x,0,t(i)));
end
%show value to console
err = expected_value - y;
figure;
subplot(2,1,1);
plot(t,y,t,expected_value);
% axis([0 1 -2 2]);
title('Signal with Trapezoidal Integration');
xlabel('time (seconds)');
legend('Observed','Expected');
subplot(2,1,2);
plot(t,err);
title('Error with Trapezoidal Integration');
xlabel('time (seconds)');

%compute continuous matlab integration
% syms x;
% fx = 2*pi*f*cos(2*pi*f*x) ;
% % fx = 2*pi*f*cos(2*pi*f*x) + 2*pi*f1*cos(2*pi*f1*x) + 2*pi*f2*cos(2*pi*f2*x)...
% %     + 2*pi*f3*cos(2*pi*f3*x) + 2*pi*f4*cos(2*pi*f4*x) + 2*pi*f5*cos(2*pi*f5*x)...
% %     + 2*pi*f6*cos(2*pi*f6*x) + 2*pi*f7*cos(2*pi*f7*x) + 2*pi*f8*cos(2*pi*f8*x)...
% %     + 2*pi*f9*cos(2*pi*f9*x) + 2*pi*f10*cos(2*pi*f10*x);
% expected_value = double(int(fx,x,0,t_final))
end